/* ========================================================================
 * sim_interface.c
 * The simulator interface for the implementation of 
 *  generalized state space Markov chain models.
 * 
 * Andr� Pedro (c) apedro@gmx.com
 * 
 * Discrete event system simulator based on a Markov chain model
 *   Simulates a generelized semi-Markov process
 *
 * ======================================================================*/
/* $Revision: 1 $ */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include "mex.h"
#include "sim_mc.h"
#include "dlinkedlist.h"

#ifndef D_LINKED_LIST
#define D_LINKED_LIST
GEN_DLINKEDLIST_FUNCTIONS(dll_sim, DLList, DLNode, double);
#endif


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, 
  const mxArray *prhs[]) {
    
    int i;
    
    int n_S, n_E, r_ES, c_ES, r_p, c_p;
    
    mxArray const *SData, *EData, *ESData, *pData, *GData, *nsimData, *seedData;
    double *S, *E, *ES;
    double *p, seed;
    int nsim; /* Number of simulations */
    
    mxArray *CData, *pathData;
    unsigned int *path;
    double *C;
    
    int n_evclk, _evclk, n_par, _par;
    double *dparameter;
    mxArray *nameev_cell, *parev_cell, *par_cell;
    char *distname;
    dparam *dp;
    
    time_t curtime;

	/* Get the current time.  */
	curtime = time (NULL);
  
    printf("GSMP simulator. Started at %s", asctime(localtime (&curtime)));
    printf("Simulation parameters:\n");
    
    if (nrhs < 7)
        mexErrMsgTxt("The number of arguments are not correct.");
    
    /* S array - States id */
    SData = prhs[0];
    S = mxGetPr(SData);
    n_S = (int)mxGetN(SData);
    printf("\t n_S: %d\n",n_S);
    
    /* E array - Events id */
    EData = prhs[1];
    E = mxGetPr(EData);
    n_E = (int)mxGetN(EData);
    printf("\t n_E: %d\n",n_E);
    
    /* ES array - for all s: for all i: Event(i) in state(s) */
    ESData = prhs[2]; /* Matrix */
    ES = mxGetPr(ESData);
    r_ES = (int)mxGetM(ESData);
    c_ES = (int)mxGetN(ESData);
    printf("\t r_ES,c_ES: %d %d\n", r_ES, c_ES);
    
    if (r_ES!=n_S)
        mexErrMsgTxt("The size of ES matrix not match.");
    
    if (c_ES!=n_E)
        mexErrMsgTxt("The size of ES matrix not match.");
    
    /* Probabilistic matrix - P(s' | s and e*) */
    pData = prhs[3];
    p = mxGetPr(pData);
    /*for (z=0; z<n_E; z++)
        for (j=0; j<n_S; j++)
            for (i=0; i<n_S; i++)
        mexPrintf("T%d: %f \n", i, p[(i+(j*n_S))+(n_S*n_S*z)]);*/
    
    r_p = (int)mxGetN(pData);
    c_p = (int)mxGetM(pData);
    printf("\t r_p,c_p: %d %d\n", r_p, c_p);
    
    /* Clock distribution matrix */
    GData = prhs[4];
    n_evclk = (int)mxGetM(GData);
    printf("\t r_G,c_G: %d %d\n", mxGetN(GData), n_evclk);
    
    DPARAM_CREATE(dparam, dp, i, n_evclk);
    
    if (n_evclk != n_E)
        mexErrMsgTxt("The number of clocks not match with the number of events.");
    
    for (_evclk = 0; _evclk < n_evclk; _evclk++) {
        /* Event distribution name */
        nameev_cell = mxGetCell(GData, _evclk);
        if (mxIsChar(nameev_cell) != 1)
            mexErrMsgTxt("Input cell must be a string.");
        distname = mxArrayToString(nameev_cell);
        
        if ((strlen(distname)+1)>5)
            mexErrMsgTxt("Distribution id much bigger.");
        sprintf((char *)dp[_evclk].nameid, "%s", distname);
        
        /* Event distribution  parameters values */
        parev_cell = mxGetCell(GData, _evclk + n_E);
        n_par = (int)mxGetN(parev_cell);
        
        for (_par=0; _par < n_par; _par++) {
            par_cell = mxGetCell(parev_cell, _par);
            /* Distribution parameters allocation in 'dp' array of structures */
            dparameter = mxGetPr(par_cell);
            dll_sim_insertbeginning(dp[_evclk].param, dll_sim_newnode(*dparameter));
        }
    }
    
    nsimData = prhs[5];
    nsim = (int) mxGetScalar(nsimData);
    printf("Number of simulations: %d\n", nsim);
    
    seedData = prhs[6];
    seed = (double) mxGetScalar(seedData);
    
    CData = mxCreateDoubleMatrix(nsim, n_E+1, mxREAL);
    C = mxGetPr(CData);
    pathData = mxCreateNumericMatrix(nsim, 2, mxUINT32_CLASS, mxREAL);
    path = (unsigned int *) mxGetData(pathData);
    
    /* Simulator */
    sim_mc(nsim,seed,S,n_S,E,n_E,ES,p,dp,C,path);
    
    printf("GSMP simulator. Terminated at %s", asctime(localtime (&curtime)));
    
    /* Outputs of mex function */
    nlhs=2;
    plhs[0]=CData;
    plhs[1]=pathData;
}
